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ABSTRACT 


The calculated stresses and displacements induced in anisotropic 
plates by short duration impact forces are presented in this report. 

The theoretical model attempts to model the response of fiber composite 
turbine fan blades to impact by foreign objects such as stones and hail 
stones. In this model the determination of the impact force uses the 
Hertz impact theory. The plate response treats the laminated blade as 
an equivalent anisotropic material using a form of Mindlin's theory 
for crystal plates. The analysis makes use of a computational tool 
called the "fast Fourier transform". Results are presented in the 
form of stress contour plots in the plane of the plate for various 
times after impact. Examination of the maximum stresses due to im- 
pact versus ply layup angle reveals that the ±15° layup angle gives 
lower flexural stresses than 0°, ±30° and ±45° cases, for 55% graphite 
fiber-epoxy matrix composite plates. 
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SUMMARY 


This report summarizes the recent progress in the attempt to model 
foreign object impact of fiber composite fan blades by small objects 
such as stones and hardstones. The high speed impact of these objects 
with composite materials should result in short duration impact times 
(< 100 ysec) . In the present model the composite blade is 
represented by an anisotropic plate of infinite extent. Thus only the 
mechanics in the area of the impact point are studied. The effect of 
edges and stress wave reflection are not treated in this report, though 
work is proceeding on this aspect of the problem. 

The combined impact contact dynamics and plate response are separated 
into two sub-problems. The. impact force time history is determined from 
the Hertz contact theory for a half space, while the plate response is 
determined from an assumed impact pressure distribution over the impact 
contact area. Linear elastodynamics is used, neglecting viscoelastic- 
plastic, nonlinear and fracture effects. Thus it is hoped that the model 
will predict, prefracture or predamage stresses. 

Using an approximate theory of anisotropic plates due to Mindlin, 
five waves are shown to make up the main part of the motion. Two of the 
waves involve inplane displacements while the other three waves involve 
flexural plate motion. The analysis makes use of a computational tool 
called the fast Fourier transform. Several computer programs were 
developed to calculate the stress levels behind the wave fronts due to 
a specified impact force distribution. The output of these programs are 



in the form of stress contour plots. It was generally found that the 
stress levels were highest in directions along the fibers. Also the 
maximum stresses appear to propagate in the lowest flexural mode. 

These waves are found to be highly dispersive and change their shape 
as the wave propagates. Examination of the maximum stresses due to 
impact versus ply layup angle, reveals that the ±15° layup angle gives 
lower flexural stresses than 0°, 30°, and 45° cases. The flexural 

stresses for the 15° case are 34% lower than those for the 45° layup 
angle case. For the interlaminar shear stresses, the values seem to 
be insensitive to layup angle. For the average membrane stresses due 
to inplane motion, the values immediately after impact appear to be 
lower for the lower fiber layup angles. The flexural stresses are 
found to depend on the ratio of impact object radius to plate thickness. 

Continued work is in progress on the edge impact problem and the 
effect of leading edge protection. 
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INTRODUCTION 


This research has been motivated by recent attempts to study the 
impact resistance of fiber composite materials. These materials are 
being considered for application to jet engine fan or compressor 
blades and must withstand the forces of impact due to the ingestion 
of objects such as birds, stones, or hailstones at speeds up to 
500 meters per second. Recognizing, of course, that inelastic 
deformation will occur at the impact point, we none the less pursue 
an elastic analysis as a prelude to the more difficult task of 
predicting permanent damage. In this report we are interested 
principally in how the energy propagates away from the impacted 
area. It has been observed that damage in these fiber composite plates 
occurs away from the impact area near edges and boundaries as well as 
at the impact point. 

In a previous paper (ref. 1) , the author presented a mathematical 
model for stress wave propagation in anisotropic plates based on the 

work of Mindlin and co-workers. Five partial differential equations 

r 

of motion were obtained for orthotropic symmetry in which the inplane 
and flexural motion were described. The two-dimensional velocity and 
wave surfaces were presented and the principal vibratory direction 
of particle motion for each wave normal was presented. 

Section II will present an analysis of the two-dimensional 
impact problem. This analysis is based on the use of a Laplace 
•transform on time and a two-dimensional Fourier transform on the 
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space variables. The solution permits the analytical inversion of 
the Laplace transform while a computational tool called the "Fast 
Fourier Transform" (ref. 2) is used to numerically invert the Fourier 
transform solution. Estimates of the impact time and force are made 
using the Hertz impact theory. 

As a preliminary to the two-dimensional problem, the author tried 
this analytical-computational technique on a few simpler, one-dimensional 
wave problems. This section will present the results of that study. 
Displacement and stresses are calculated for a short duration line 
force normal to an anisotropic plate. The responses for various 
fiber layup angles are compared. 
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SECTION I: ONE-DIMENSIONAL TRANSIENT WAVES IN ANISOTROPIC PLATES 


The mathematical model used in this paper is called an effective 
modulus theory of composites. The model is valid so long as the scale 
of the changes in stress levels (e.g. wavelength) is much larger than 
the sizes of the composite constituents (e.g. fiber diameter, ply 
spacing, etc.). This assumption has been used by many researchers 
to derive equivalent elastic constants from long wavelength or long 
pulse length wave propagation tests (e.g. Tauchert and Moon, (ref. 3)). 
Equivalent moduli determined in this way have checked very closely with 
elastic constants obtained from static and vibratory tests. Tauchert 
and Guzelsu (ref. 4) , have used ultrasonic tests to measure dispersion 
in composites and found no significant departure from the effective 
modulus theory for longitudinal waves in boron fiber/epoxy matrix 
composites until the frequency exceeded 5 megahertz. Shear waves, 
however, began to show dispersion at about 0.5 megahertz. This 
corresponds to a wavelength to fiber diameter ratio of about 40. 

Sun et al . (ref. 5), in calculation on laminated composites, reports 
similar results. Thus the effective modulus theory used in this paper 
is subject to the restriction that the impact pulse spectrum have its 
most significant wavelengths greater than about 100 times the fiber 
diameter. Further discussion of dispersion in composites may be found 
in a review of the subject by the author, (ref. 6). 

Equations of Motion 

In the approximate theory of elastic plates in this study (ref. 7) , 
the displacement is expanded in the thickness variable by using Legendre 
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polynomials P (n) , 


u = E P (n)u 
^ n=o n ^ 


(n) 


(x ,x ,t) 
1 3 


(I-D 


where n = x /b ,b is the half thickness, (see Fig. 1-1) . A variational 

2 

method is used to obtain the equations of motion. For the purpose of 

this work, the series was truncated at n = 2, with only seven terms 

C2'\ 

used; namely, u° , u° , u°, u 1 , u 1 , u 1 , u L J . This leaves all the strains, 
1 2 3 1 2 3 2 

except e , e , as linear functions of the thickness variable n . 

12 32 

To further simplify the equations, the inertias and high frequency 

terms (derivatives) of u 1 , u^ were dropped, resulting in explicit 

2 2 

expressions for these terms ; 
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where (-q ) is the transverse loading on the top plate surface. Elimination 
2 

r 2^ 

of the terms u 1 and u^ 

2 2 

results in the equations 


from the remaining five differential equations 
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In his approximate theory of crystal plates (ref. 7), Mindlin uses 
correction constants to adjust the thickness vibration mode to match the 

exact theory. Thus he replaces C and C by' k C and k C respectively 

44 66 1 44 3 66 

In our case C = C and k = k = tt 2 /12 . 

44 66 13 

The stresses in the plate for orthotropic symmetry take the following 

form: 
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The interlaminar shear stresses t and t are quadratic functions 

12 32 

of x and cannot be accurately found from u . However in analogy 

2 ^ 

to classical plate theories, we can integrate the original stress 
equations of motion to obtain 


t 
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The in-plane stresses t , t , t are comprised of two parts. 

11 33 13 

One part is uniform through the thickness and is governed by dispersion 
free equations (I-3a) and (I-3b).. The flexural part of the stress is 
linear in the thickness variable and is governed by equations (I-4a) 
through (I-4c) which exhibit strong dispersion of wave phenomena. 
One-Dimensional Wave Solutions 

We imagine an infinite plate on which a distributed load is placed 
along the line jj • £ = constant (see Fig. 1). The vector j\ is 
normal to the line, and we assumed the load on the plate is independent 
of the position along the line. Thus the surface load has the form 

q(x ,x ,t) = q Oj - £>t) 

Because of this relation, the motion is also assumed to have this form, 
i . e. 

u°(x , x ,t) = u°(n.r,t), etc. 

113 1 ^ ^ 

For derivatives of functions of this form, the following expressions 
are useful ; 

9f 0 ' 9f 0j-£> * 

— - = f'cos a, — = f'sin a 

9x 9x 

1 3 
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where primes indicate differentiation with respect to the variable 


5 = S * £ » and 




x cos a + x sin a 
1 3 


The angle a , denotes the departure of j\ from the symmetry 

axis x = 0 . To find solutions to the equations (1-3) and (1-4), 
3 

one writes the unknowns in the forms 


u ° = u ° = 


(1-7) 


u 0 = WOj.^.t), uj = - b ^(jv^t), u* = b 'i' i (?J.^,t) 


Note that q is a prescribed function in time, along a line parallel 
to the normal n . This reduces the differential equations of motion 
to equations in two variables £ , t where £ = n • r . The formal 
solution to this initial value problem is found by taking a Laplace 
transform on the time variable and a Fourier transform on the space 
variable. These operations are defined as follows: 


L[f]=f(s)=J e _St f(t)dt 
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T [f ] = f (k) 


/ e lk f(Od? 


( 1 - 8 ) 
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Thus when the displacements have the forms given by (1-7) and 
the operations (1-8) are performed on the equations of motion, the 
following algebraic set of equations result, where initial values 
of }i are assumed to be zero and where q = 0 for t < 0 . 
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where 


= C cos 2 a + C sin 2 a 
11 11 55 


A - C sin 2 a + C cos 2 a 
22 33 55 

A = C + C 
12 55 13 

= 3 C /b 2 p 
0 66 ■ 


0 is the frequency of the first thickness shear mode from the exact 
theory. Also, 


a = C /2C 
1 12 22 


a = C / 2C 
2 23 22 


a = C /2bC , a = C /2bC 
3 12 22 4 23 22 


Inversion of Laplace Transform 

When q(£,t) is given, these equations may be solved in the 

a, 

transformed variables. For example consider the following case. 

CASE I.) Midplane Motion; Normal Load (q = q =0) 

1 3 

^ A 

U = iq {(A +s 2 /k 2 )a - A a sin 2 a }cosa/kA 

1 '2 22 1 12 2 


0 - iq ((A +s 2 /k 2 )a - A a cos 2 a}sina/kA 

3 2 11 2 12 1 


where 


( 1 - 11 ) 


A - (A + s 2 /k 2 ) (A + s 2 /k 2 ) - A 2 cos 2 asin 2 a 
11 22 12 

These solutions in the transformed plane are thus of the form 
0 = iq P(s 2 )/kA(s 2 ) 


12 



where P(s 2 ), A(s 2 ) are polynomials in s 2 . A(s 2 ) has four zeros in the 

complex s plane for each k. The two roots |s /k| and |s /k| 

l 2 

correspond to the two real wave speeds associated with the propagation 
of plane waves in the plate (ref. 1). In fact, , is simply 

the two-dimensional acoustic tensor for the plate as discussed in a 
previous paper. We are thus guaranteed four pure imaginary roots 
for s 


s = ±iv k , ±iv k 

i 2 


The inversion can then be done by use of the convolution theorem; 


i.e., if 


G(t) = L - 1 [ ] 

A(S 2 ) 


U(k,t) = L _ 1 [U] = f t q (k ,t-x) G(x)dx 

r\ 2 


For our case it can easily be shown that 


G(t) = - 


kP(-k 2 v 2 ) sinkVjt kP(-k 2 v 2 ) sinkv 2 t 


(v 2 -v 2 ) 
1 2 


(v 2 -v 2 ) 
1 2 


CASE II.) Flexural Motion; Normal Load. 


In a similar fashion when q is given, W, f or ? can be 

1 3 

solved and have the form 
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However, in this case the roots of A = 0 are not proportional to k . 
This means that the velocity of wave propagation depends on the wavelength. 
As k -* °° it is known that the wave speeds become independent of k and, 
in fact, two of them equal the non-dispersive wave speeds v , v found 
in the previous example (ref. 1) . It is also known from the dispersion 
relations for plates (ref. 7), that for k real, there will be six 
pure imaginary roots of A (s 2 ) = 0 at 

s = ±ia) (k) , ±i<jj (k) , ±iw (k) 
l 2 3 

The relations w^(k) are known as the branches of the dispersion relation 
for these plates. One branch goes through the origin, i.e. 

w (k) ->0 as k -*■ 0 

1 

The other two branches are sometimes called optical since as the 

wavelength becomes infinite, i.e. k -*■ 0, to and to denote the vibra- 

2 3 

tory frequency of the thickness shear mode for the plate 

(d (0) = u> (0) = 0 
2 3 0 

A typical dispersion curve is shown in Figure 2. 

The inversion again makes use of the convolution theorem and the 
residue calculus to invert 1/A (s 2 ). Thus we obtain 



H(t) = 


sincjjt 


+ 


sinw 2 t 


(1-15) 


RO*^) 
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2 1 3 1 
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12 3 2 


0) 
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+ : — — — - 

(w 2 -to 2 ) (u ) 2 -(jO 2 ) 

13 2 3 


sinu 3 t 


0) 
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Numerical Inversion of the Fourier Transform 

It is at this juncture that the solutions to most problems in 
transient wave propagation are limited by analytic skill in extracting 
information from the inversion. In a few cases the complete solution 
can be given, while in most, far field, near field, short and long 
time approximations must be invoked. The numerical inversion of the 
transform has until recent years been limited by the calculation of an 
integral for each point in the space £ . However, recently a technique 
has been developed which obviates the need for a separate inversion 
for each point in the real space ? . Known as the "Fast Fourier 
Transform" (ref. 8) , it takes a sampling matrix of the transform of U, 

A A A 

say (U(k^), U(k ), ... Ufk^)) and returns a sampling matrix of the 
original function 


(UU ), U(c ), U(c ), ... UU N )) 
12 3 N 


For a large enough . N, one will obtain a picture of the original function 
U(c). Application of this algorithm to our problem is described below. 



The specific operation that these computational algorithms perform 
is called a finite, discrete Fourier transform. Thus if D(I) is a 
one x N matrix of data, the output of the algorithm T(J) is given by 


T (J) 


1 = 1 N 


(1-16) 


The properties of this operator are similar to the conventional Fourier 
operators in the space of continuous variables (ref. 8). The sum is 
not performed in a direct manner on the computer but in a manner which 
reduces the number of arithmetic steps and makes the operation attractive 
as a computational tool. 

In the solution of partial differential equations by Fourier 
transforms we are led to evaluate integrals of the type 

U( C) = j T~ U(k)e" lkC dk 

^ -00 

If significant changes in U(c) occur over distances greater than 
A , then the largest wave number of interest will be 

k = 2t:/X 

Thus we may be satisfied with an approximation to U(C) of the form 

fiU) = 2 I* U(k)e' lkc dk 
Zl T 

-K 
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or shifting the coordinates 


'V 

U = 


2tt 


! 2k U(k-K)e" lkC dk 


This latter integral may further be approximated by the limit of a sum 


ft = 


K iK? 
ttN e 


E 


U(k I -K)e“ lk I ? 


(1-17) 


where 


or 


k = — 

I N 


+ (I-D 


4 TT 
NX 


I + f 1 - 1 ^ 


ft = e 1 K ^^ 1 - N' ) Z U(kj-K) e' 12 ^ 1 " 1 ^ ^ 


So far U(c) has remained a continuous function of C . However at 
the points 


Cj = j(J-l) , J = 1, 2, .. N 


the summation becomes identical to the finite descrete Fourier transform 
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T (J) 


(1-18) 


&(Cj) = ?I(J) 


iir (J-l) (N-l) k 
6 N irN 


where T(J) is defined in (1-16) and 
D C I ) = U(k T -K) 


Note that this scheme gives details of U(£) of distance no smaller than 
A/2 . This was implicit in choosing maximum wave number. Further when 
N is even, the output matrix U(J) describes the function only up to 
S = - 1) . When is real, and N even, we have the 

identities 

ft (j + 1) = 0 , Imag (U C J) ) = 0 

N <v N 

U(j + 1 + <) = - U(j + 1 - K) 

Thus the one-dimensional output matrix b(j) which approximates our 
original function is antisymmetric about J = y + 1 and symmetric 
about J = 1 . Since lj(J) = 0 for J = 1 + N/2 we must choose 

N large enough, or the time small enough to ensure that our wave has 
not reached £ = A(N + 2)/4 . In this approximation, we have 
effectively replaced a single impact source by 
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an infinite periodic array of sources and negative sources. Our 
approximation will be valid as long as the waves from each of these 
sources do not interact with each other. 

Notes on the Numerical Fourier Inversion 

There were several checks made of the fast Fourier computer 
routine (ref. 8) used in this paper. First, known functions were 
transformed analytically and inverted numerically. The results of 
these tests revealed that the inversion program works best on continuous 
functions if one is to avoid spurious oscillations near points of 
discontinuity . 

Second, the initial value problem for a nondispersive medium 
such as a string was programmed. The results of this test are shown 
in Figure 3a. The shape of the string at the initial time and subsequent 
times as predicted by the computer preserve the original triangular 
shape and exhibit the wave shift at the proper speed. 

Lastly the impact of a string was programmed using the same force 
applied to the plate problem. The theoretical result predicts a constant 
displacement behind the wave front as shown in Figure 3b. 

For the line impact of an anisotropic Mindlin plate, a specific 
load distribution was chosen for ease of analytical calculation of the 
Fourier transform and Laplace convolution integral; 


for 


for 


u < a and t < x 

o 

n 1 n TT^-i . Wt 

q = -P -=-(1 + cos—) sm — 
? o 2 a x 

I cl > a or t > x 


q 2 = 


( 1 - 19 ) 


0 
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By using the fast Fourier transform (with a proper sign change) 
one could find the transform for an arbitrary load distribution. 

However, this has not been done. Nor has the author used the load 
distribution for a Hertz contact problem, as might be supposed in an 
impact problem (reasons for this will be given below). However the 
chosen distribution, it is hoped, will exhibit most of the salient 
effects of the impact of a plate. 

While ad hoc, the particular choice of the load distribution is 
not entirely arbitrary. Continuity of the first and second derivatives 
(as the function (1-19) exhibits) is dictated by a desire to have all 
stresses continuous (i.e. to avoid shocks) and thus spurious oscillations 
in the numerical inversion. This stems from the fact that in the 
Mindlin theory the midplane stresses, having as wave sources the 
term Vq 2 , have the following form for their transforms 

~ rt 

t 0 'v, J kq (k, t) sinkv(t-x)dx , (a, 3 = 1,3) 

“P o 2 

where v denotes either the first or second wave speed. Thus when 

A 

q 2 has the form 

q = Q(k)S(x) 

2 


the stresses are proportional to 


t 'v k Q(k) sinkvt ^ T 


where. 


9 2 O' 

h 2 j 


T [F(vt-c)] 


F(vt-c) =1, U | < vt 
= 0 Id < vt 
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Thus for continuous stresses at early times after the wave arrival, the 
spacial part of must have continuous second derivatives, which led 

to the choice of (1-19) . This conclusion was found also in the numerical 
results when non-smooth load distributions were used. 

Results for the Line Impact Problem 

Using the transient load distribution described above, calculations 
were made, on an IBM 360-91 computer, of the stresses and dis- 
placements in both a classical and Mindlin plate due to impact type 
loading. The results are shown in Figures 4-8. 

One of the important parameters in this problem is the ratio of 
load extent to plate half thickness, a/b. When a/b is of order 
unity or less, one would expect that shear effects would become important 
and the Mindlin and classical plate solutions would differ. This is 
so ,as shown in Figure 4 for the case of a/b = 1 . However in Figure 5 
when a/b = 10 the displacements calculated by both Mindlin and classical 
theories do not differ by very much. 

The remaining figures are for the Mindlin plate and were calculated 
for the composite, 55% graphite fiber-epoxy matrix,' using equivalent 
elastic constants for various layup angles of the fibers. The constants 
were taken from Chamis (ref. 9) . 

In Figure 6 we have plotted the center plate displacement versus 
time for both the classical and Mindlin theories. For the classical 
plate this function can be found explicitly. The displacement increases 
as the square root of time when the impact force is a delta function 



in time and space. This result is confirmed by the numerical results 
in Figure 6 and is also the case for the Mindlin plate for large time. 
(Note that for a string, the displacement at the center is time 
independent after impact.) 

The effect of layup angle on the plate response is shown in 
Figure 7, for various times after impact. As can be seen, the 
displacement is somewhat insensitive to layup angles of up to 

about ±15° for either line loads along the x axis or along the 

3 

,x axis . 

1 

Finally in Figure 8 the induced membrane or average in-plane 

stress 1/2 (t + t ) is shown for various layup angles at a 
11 33 

time immediately after the impact time. One can see that the 

initial compression pulses are preceeded by wave fronts which 

vary in magnitude and velocity with layup angle. Also for a 

wave propagating along the x axis (load on the x axis) , 

1 3 

the ±45° layup results in higher stresses than for other layup 
angles. While the flexural stresses are higher than the membrane 
stresses, the membrane stresses will propagate ahead of the bending 
waves and will have a tension pulse in the signal which might 
cause splitting through the plate or plys as have been observed 
in experiments . 

In the case of one-dimensional wave propagation in a direction 
n = [cosot, sina],, the direction of the displacement is not parallel 
to ji but is known for the particular wave in question. The direction 
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of displacement for the first two in-plane waves was given in Figure 6 
of Ref. 1 

[u°, u°] = U(n • r-vt) [cos 8, sin 3] 

13 % 


From (1-5), the stress may be determined. In particular, the average 

in-plane stresses t° , t° , t° can be found as functions of the 

11 33 13 

mean stress a = l/2(t° + t° ) ; 

1 1 33 


t° = 2a(C cos a cos 3 + C sin a sin 8)/D 

11 n 13 


t° = 2a (C sin a sin 8 + C cos a cos g)/D 


3 3 33 


1 3 


tP = 2 a C (cos 8 sin a + cos a sin g)/D 
1 3 55 


where 


(C +C ) cos a cos 8 + (C + C ) sin a sin 8 
11 13 33 13 
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SECTION II : STRESS WAVES DUE TO CENTRAL IMPACT OF ANISOTROPIC PLATES 


In Section I of this report, an analytical- computational 
method was described for the calculation of the impact induced stresses 
behind the wave front. The method was checked for known, one-dimensional 
problems such as the impact of a string, and the line impact of a classical 
anistotropic plate. The induced stresses were given in terms 
of the impact pressure. However, no prediction was made of the maximum 
impact pressure in terms of the velocity and mass of the impacting body. 

In this part of the report, the problem of normal or central impact 
is discussed. Two-dimensional stress patterns are presented in terms of 
the impact pressure and an estimate is made of its magnitude for a 
spherical impacting body of known velocity and mass. 

Description of the Mathematical Model 

Before discussing the results, a discussion of the assumptions made 
in the model will be given. In practice, fiber composite plates are made 
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up of a number of unidirectional plys oriented at various angles to 
obtain certain desired properties. When the properties and lay up 
angles of each ply are known, one can obtain the plate constants by 
averaging across the thickness. These average elastic constants 
represent an equivalent anisotropic-homogeneous material. In the present 
work the laminated plate is replaced by an equivalent anisotropic plate. 

The equivalent elastic constants were obtained from a computer code 
developed by Chamis (ref. 9) . The following additional assumptions are 
implicit in this model. 

Linear Elastic Properties 

This model is based on linear elasticity. Thus, plastic, viscoelastic, 

and fracture effects are not taken into account. The results of these 

effects is to decrease the amount of wave energy that can propagate into 

the plate. Thus, the elastic analysis represents an upper bound 

on the actual stresses during impact. The effect of initial stress in 
the plate has also not been considered in this report. 

Boundaries and Tapering 

Another limitation of this work is the neglect of the effect of 
boundaries or edges of the composite fan blade. The reflections from 
boundaries can of course be handled with the linear elastic theory. 

However, if one is interested in the maximum stresses, these should occur 
near or at the impact point except for the case of edge or near edge 
impact. Also, for small impact times, the waves behave as if the finite 
plate were infinite. The case of edge impact is under study at this 
writing. 

In this report the plate thickness is assumed to be constant, while 
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I 


in a fan blade taper and twist will be present. The neglect of 
these effects seem to be of lesser importance than those due to 
inelasticity and boundaries. 

Thickness Reflections and the Mindlin Theory 

The use of an equivalent anisotropic plate neglects the reflections 
of waves at the ply boundaries. However, this approximation will be valid 
for wavelengths* greater than the ply thickness, and valid also for 
wavelengths greater than the plate thickness. There are analyses which 
consider wave-ply interactions, (Ref. 10), but in general the change 
in stress must occur over distances comparable to ply or fiber thickness 
for these effects to become important. 

In the Mindlin theory of plates , the wave reflections are averaged 
through the thickness, (ref. 7). In dynamic loading of a plate, the stresses 
will propagate through the plate thickness as well as away from the 
impact point. These waves will suffer many reflections as they propa- 
gate back and forth between the plate surfaces. The calculations of these 
many reflections for a three-dimensional plate thus become impos- 
sible for anything but very early times after impact. The Mindlin plate 
theory thus restricts the mathematical equations to a description of 
the average plate midplane motion and rotations. 

In the Mindlin theory used in this report the plate displacements 
u = (u^, u^y u^,) (See Fig. 9) have the form (See Ref. 1) 


*An effective wavelength can be defined as the contact time of impact 
multiplied by the smallest wave velocity in the material. 
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Cii-l) 


u = u°(x 1> x 3 ,t) + Jp u^ (x^jX^ ,t) + y[3(^-) 2 -l] u^ (x^x^t) 


The equations governing the motion are listed below. The first two 

O O 

l l’ U 3 


equations are for the inplane displacements u° ,u° and the next three 


govern the transverse displacement u° and flexural rotations u^/b , 
u 3 /b . The transverse impact loading is represented by 
which we will prescribe in the following section. 
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where 
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12 

22 

2 

32 
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C 

1 3 


C 

13 




kC 


66 


The factor k is a correction factor introduced by Mindlin and 
equal to k = tt 2 / 12 

If the plate is struck at a point, the energy will propagate into 
the composite in the form of elastic waves (see Ref. 1). Two of these 
waves will have anisotropic wave fronts and will depend on the fiber 
layup angle of the composite as shown in Figure 10. The average in-plane 
stresses across the thickness will propagate in this manner. The 
flexural energy will propagate behind a slower isotropic wave front. 

The stresses in the plate are given in terms of the displacement 
(Section I of this report) . The mathematical problem to be solved 
consists of the following: find a solution to the coupled partial 

differential equations (II-2) , (II-3) when the loading function q is a 

2 

prescribed function of the plate coordinates x ^> x 3 and time, and the 
plate is initially at rest. 

The solution of this problem was accomplished using a combination 
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analytical and computational technique involving Laplace and Fourier 
transforms. Details of finding the solution are discussed in Appendix A 
of this report. Results of these calculations are discussed below. 

Before considering the results, we must discuss the choice of the 
impact loading function q 2 • 

To determine the maximum contact pressure one must know the 
impact pressure distribution on the plate surface. In the static Hertz 
theory, this pressure is given by (see Ref. 11 and 12) 

q = -P n (l " (f ) 2 ) 1/2 C 11 - 4 ) 

2 U a 

This distribution is not suitable for our theory since the infinite 
slope of P(r), at r = a , would appear in our approximate equations 
(II-2) and (II-3) as a source term. In fact, it is shown in Section I of this 
report that the second derivative of P(r) must be finite at r = a 
for the stresses to be continuous. To resolve this problem, an ad hoc 
pressure distribution is assumed which produces finite stresses; 


q = -P (1 - 2(f) 2 + (f) 4 )f(t); r < a 
2 0 a a — 

q = 0 , r > a 
2 


C 1 1 - 5 ) 


f(t) = sin 7Tt/x , t <_ T , f(t) =0 t > T 
It should be noted that in an actual impact, the impact area 
changes with time. This effect is difficult to model in the analytical 
part of the problem. Instead we have chosen to prescribe a time varying 
pressure over the maximum contact area . 
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The total force produced by such a pressure distribution is 


F^= 2 tt\ q(r) r dr 


'0 

7ra 2 P 
F = ^ 

0 3 


Cl 1-6) 


Using ( 9 ) we may calculate the total impulse to be 


, T 


I = \ F sin— dt = ~r a 2 P x 

0 x 3 0 

0 


CH-7) 


If the rebound velocity is V , and the inital velocity V , 

1 0 

the impulse given by 

I = M(V +V. ) C 1 1 - 8) 

10 

Herztian Impact Stress 

In order to determine the impact stress equations (II-2, 3), the 
impact pressure distribution must be found in terms of the object mass, 
velocity, geometry and elastic properties. No exact solutions to this 
dynamic problem are known in the theory of elasticity. However the 
static theory of contact of elastic isotropic bodies is known and was 
developed by Hertz (ref. 11). It is known, however, that even under 
small contact forces the contact pressure exceeds the elastic limit 
of conventional materials (see e.g. Goldsmith, (ref. 12)). One would 
expect then that for high speed impact (>100 m/sec.) the local contact 
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stresses would exceed the elastic limit of most materials. In spite of 
this violation of the elastic assumption of the Hertz theory of 
impact, experiments have shown that the impact times predicted 
by the elastic theory agree reasonably well with data from impact 
tests. The stresses predicted by this theory however are upper 
bounds on the actual stresses. Also if the surface of impact belongs 
to a structure which can move, then the contact times predicted by 
this theory are lower bounds on the actual time of contact of object 
and structural surface. With these restrictions in mind the Hertz 
theory of impact will be reviewed as it can be adapted to composite 
materials . 

Consider the contact of a sphere of radius R and elastic half 
space. The contact force between the two bodies F, is related to 
the relative approach of the bodies a, and has the dependence 




(II-9) 


where is a constant dependent on the properties of the bodies. 

When both bodies are isotropic this constant is given by 


k 

2 


4 R l/2 
3 R 




(II- 10) 


where v, E are Poisson's ratio and Young ' s modulus respectively for 
each body. 

Composite materials however are anisotropic in general. The 
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corresponding problem for the contact of anisotropic bodies has 
recently been given by Willis (ref. 13). For this case the contact 
area is an ellipse, the dimensions of which must be determined from 
algebraic equations. The shape of this ellipse for typical composite 
anisotropy has not been determined, though from experiments of the 
author on graphite/epoxy and boron/ aluminum composites seem to show 
that the contact area is close to a circle. It is assumed then that 
for the impact of an isotropic sphere with a composite surface the 
contact area is a circle. Also the constant k is replaced by 


k 

2 


iR 1/2 


!-vf 



( 11 - 11 ) 


where E ,v are the elastic constants of the sphere, and C is the 
1 1 22 

elastic modulus of the composite plate. This assumption of course is 
just an educated guess. 

In the Hertz impact theory the force (II-9), which was determined 
from a static solution, is used in Newtons law for the sphere 




( 11 - 12 ) 


where M is the mass of the sphere and V the instantaneous velocity of 
the sphere. 

This assumption is only valid when the contact time is much greater 
than the time for elastic waves in the sphere to traverse the object. 



Further when the motion of the plate becomes large during the contact 
interval, (11-12) must be replaced by 


M — = k(U-W) 3/2 (11-13) 

dt 2 


where U is the displacement of the sphere and W is the displacement 
of the plate at the point of contact. This problem requires simultaneous 
determination of the motion of the sphere and plate. This has not been 
done in this report but work on this problem is in progress. 

When W — 0 during the contact period, the time of contact can be 
shown to be (see e.g., Goldsmith (ref. 12)) 


T 



da 

[V 2 - ^ « 5/2 1 I/2 


2.943 Ql 
_ 


(11-14) 


where a is the maximum approach of the impacting object. This value 
is given by 



5 

4 


2/5 


MV 2 

k 


2 


(11-15) 


For elastic impact the maximum radius of contact is given by 


a 



(11-16) 


The impact force in this theory as determined from the solution of 
(11-12) is given by 
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F = F sin irt/x , t < x 
0 ~ 

= 0 , t > x 


(H-17) 


where F =3.36 — . The pressure distribution in this theory is given • 
by tiie expression 

r „ 1/2 

P = P (1 - (f) 2 ) 

0 a 

7 

where P = — F /ua 2 . It is interesting to note that the impact 

0 2 0 

pressure P , is independent of the radius of the impacting object. 

0 

This distribution, however, is determined from a static elasticity 
solution and should not be expected to hold under dynamic impact. Of 
course the contact radius varies with time reaching a maximum value 
C I 1-4) at half the period. 

Another limitation of the theory is the fact that it predicts 

a rebound velocity equal to the initial velocity V . In an actual 

0 

impact , momentum would be transmitted to the structure, thus changing 

rebound velocity to something less than V . Solution of the coupled 

0 

problem (11-13) should enable a better determinatinn of the rebound 
velocity. 

A further improvement of the theory might be achieved if a more 
general contact law is used, say (called the Meyer Law, see (ref. 12)) 
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where n, k, would be determined from static experiments on composite 
materials. This approach is presently under investigation. 

Well aware of the long list of limitations on the above model for 
impact, we have nonetheless used the above procedure to obtain estimates 
of the impact time, contact force, and maximum contact pressure for 
different impacting bodies and speeds. The results of these calculations 
are shown in Figures 11 and 12 where data is presented for spherical 
ice particles and for granite-like stones. 

For impact speeds from 100 to 500 meters/sec. and 1/2 to 2 cm. 
diameter granite spheres, the contact times range from 15 to 85 ysec. 

In summary, these impact formuli reveal the following dependence 
on impact velocity: 


T. 'V- 1/V 


1/5 


P % 
0 


V 


2/5 


These results should only be used as guidelines, since many assumptions 

are made which break down at high velocities. 

Goldsmith and Lyman (ref. 14) have shown the Hertz 

theory to be remarkably valid insofar as contact time and peak pressure 

are concerned for the impact of hard steel spheres (1/2" diameter) into 

a hard steel surface for velocities up to 300 ft. /sec. (91.5 meters/sec.). 

The calculations in this report based on the Hertz law of impact extend 

well beyond this limit (100-500 meters/sec.). Thus the contact times 

and maximum impact pressures presented in Figures H> 12 for graph ite/epoxv 

can only be used as rough guides for which the values for contact times 

are lower bounds on the actual times and the values for pressure are 

upper bounds . 
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One further reason why the Hertz law is not valid is that for velocities 
in the range of 100 to 500 meters/sec., the contact area diameter 
approaches the diameter of the striking object ,which certainly 
violates one of the assumptions in the Hertz theory. 

At the writing of this report, work is in progress on using a 
proper Hertz contact theory for anisotropic materials. However, 
for large deformations, plates, and fracture effects in the impact 
zone, a more detailed analysis of the impact zone is needed than 
is given in this report. 

Results for Impact Stresses 

Solutions to equations (II-2, 11-3)^ which govern the central 
impact of anisotropic plates, were found for impact-like pressures 
using an analytical/computational method as described in Appendix A 
of this report. In this model there are eight different stresses (see 
Section I) associated with the membrane, flexural and interlaminar 
stresses. The presentation of all of these stresses for different 
times after impact and various layup angles becomes an enormous task. Instead 
certain key stresses or stress measures are presented in this report 
to give an overall view of the stress picture. 

The three stress measures chosen for this report were the average 
membrane stress (t° x + t° )/2 , the average flexural stress (t^ + t^)/2 
at the surface of the plate, and the maximum interlaminar shear stress 
given by (t|i + t| 3 ) 1 . These stresses are not necessarily the 
maximum stresses at a point, but they are independent of the orientation 



about an axis normal to the plate. The program also allows individual 
stresses to be obtained if desired. 

The stresses were calculated in a quarter plane of the plate for a 
specific time after the initiation of impact and were normalized with 
respect to the maximum impact pressure as calculated in the above section. 

The data is presented for various times and layup angles in the form of 

stress contour plots (Fig. 14-25) which were made on a "Calcomp" plotter 

at the NASA Lewis Research Center, (ref. 16). Superimposed on these curves are the 

theoretical wave front for the particular wave in question and the 
radius of the circle which bounds the impact pressure. The wave 
front calculations were based on the work and represent a 

locus of wave surfaces originating from the edge of the impact circle 
for a given time after start of impact. 

These results show the effect of the change of fiber layup 
angles on the stress distributions. In general, the maximum stresses 
tend to lie along the fiber directions. For most of the cases treated, 
the significant stressed region revealed by the calculations is bounded 
by the theoretical wave fronts as calculated in (ref. 1) . This provides 
a check on the accuracy of the approximations made in the numerical 
computation scheme. 
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For this elastic model, the maximum mean stresses ^(t + t ) 

I n 33 

occur at about the end of the impact time . The question about an 

optimum fiber layup angle is partially answered by the data in 

Fig. 33. For the flexural stresses , the optimum layup angle is 

about ±15° showing a 34% lower mean stress level than the ±45° 

case. However, regarding the interlaminar stresses, for the same 

impact conditions, there seems to be little difference in the maximum 

stress level with layup angle despite significant changes in stress 

distribution in space with layup angle. 

For the average membrane mean stress t° + t° , immediately 

11 33 

after impact, the lower layup angle plates yield lower maximum 
stresses. However, at later times, the ±45° layup case appears to 
give a lower maximum stress in the plate. 

Another result of these calculations is that the induced stresses 
depend on the impact circle radius- to-plate” thickness ratio (a/2b). 

On the other hand, the impact circle depends on the incoming particle 
velocity as determined by equation (8). Thus, for each impact velocity, 

a different impact radius/ thickness ratio must be chosen as well as a 
different contact time. The integration of these two programs has not 

been performed at this time but will be attempted in the near future. 

Of course, to evaluate the possibility of fracture or failure of 
the composite under impact, the complete stress matrix at a point must 
be known, as well as the failure criteria, which will itself be 
anisotropic (ref. 8) . The integration of programs described in this 
section into a computer code suitable for use for pre-design 
engineering calculations is to be the next phase of this report. 
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CONCLUSIONS 


The successful application of composite materials to jet engine 
fan blades depends, in part, on the ability of these materials to 
retain structural integrity under transient loadings due to bird 
strikes or hailstone impact or other foreign object encounter. 

While there are a number of experimental investigations connected with 
this problem, theoretical understanding of impact response and damage 
is lacking. Such understanding might enable a reduction in costly 
empirical studies and testing. This report presents the first of a 
series of analytical models to attempt to understand impact mechanics 
of composite plates. 

Using the method of stress wave analysis, the stresses induced 
during and after impact with a line load and central impact have been 
determined. The model has been put into a computer program where the 
transverse loading force on the composite plate is arbitrary. For 
central impact the results indicate that the energy propagates into 
the plate in the form of extensional and flexural waves. At several 
impact circle diameters from the center, the highest stresses propagate 
along the fiber directions. The speed of propagation varies with the 
direction in the plane of the plate. 

It has been determined that for composites similar to graphite/ 
epoxy, there is an optimum layup angle near ±15° for which the flexural 
stresses are minimized . The maximum stress in a plate without edges, 
due to central transverse impact, occurs at the end of the contact time. 
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Using a modified Hertz contact theory, an estimate of the contact 
times, pressures and forces for various impact velocities and sizes 
of ice and granite spheres has been calculated. However, the effect 
of the plate motion has been neglected, which might increase the 
calculated contact time. 

The next parameters to study in this program are the effect of 
edges on impact, and the effect of plate motion on the contact time 
and pressures. Also the validity of the Hertz theory is in question 
for the impact velocity range of interest. A modified Hertz theory 
or a fluid-solid interaction model should be developed. Some of these 
effects will be investigated in a continuing study of impact and 
stress waves in composite fan blades. 
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16 


17 


18 


19 


20 


21 


22 


23 


24 


25 


Average Membrane Stress Contours (t +t )/2 for 

11 33 0 

Graphite Fiber-Epoxy Matrix, After Impact for 15 
Fiber Layup Angle, Normalized Time t=2 T 

Average Membrane Stress Contours (t +t )/2 for 

11 33 0 

Graphite Fiber-Epoxy Matrix, After Impact for 30 
Fiber Layup Angle, Normalized Time t=2x 


Average Membrane Stress Contours (t +t )/2 for 

ll 33 0 

Graphite Fiber-Epoxy Matrix, After Impact for 45 

Fiber Layup Angle, Normalized Time t=x 

Average Membrane Stress Contours (t +t )/2 for 

11 33 o 

Graphite Fiber-Epoxy Matrix, After Impact for 45 

Fiber Layup Angle, Normalized Time t=3x 

Average Membrane Stress Contours (t +t )/2 for 

11 33 0 

Graphite Fiber-Epoxy Matrix, After Impact for 45 

Fiber Layup Angle, Normalized Time t=2x 

Average Flexural Stress Contours ft +t )/2 for 

11 33 0 

Graphite Fiber-Epoxy Matrix, After Impact for 15 

Layup Angle, Thickness/ Impact Diameter Ratio 10.0 

Average Flexural Stress Contours ft +t )/2 for 

11 33 0 

Graphite Fiber- Epoxy Matrix, After Impact for 30 

Layup Angle, Thickness/ Impact Diameter Ratio 10.0 


Average Flexural Stress Contours (t +t )/2 for 

11 33 0 

Graphite Fiber-Epoxy Matrix, After Impact for 45 
Layup Angle, Thickness/Impact Diameter Ratio 10.0 

Average Flexural Stress Contours (t +t )/2 for 

11 33 

Graphite Fiber-Epoxy Matrix, After Impact for 30 , 45° 
Layup Angles, Thickness/Impact Diameter Ratio = 10.0 


Maximum Interlaminar Shear Contours for Graphite 
Fiber-Epoxy Matrix After Impact for ±30 Fiber 
Layup Angle, Thickness/Impact Radius Ratio = 10.0 


Maximum Interlaminar Shear Contours for Graphite 
Fiber-Epoxy Matrix After Impact for ±45 Fiber 
Layup Angle, Thickness/Impact Radius Ratio = 10.0 
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26 Midplane Transverse Displacement Distribution in the Quarter 
Plane of the P^ate for Graphite Fiber-Epoxy Matrix After 
Impact for ±45 Fiber Layup Angle, Thickness/ Impact Radius 
Ratio =1.0 

27 Midplane Transverse Displacement Distribution in the Quarter 
Plane of the Plate for Graphite Fiber-Epoxy Matrix After 
Impact for ±45° Fiber Layup Angle, Thickness/Impact Radius 
Ratio =1.0 

28 Average Membrane Stress Distribution in the Quarter Plane 
of the Plate for Graphite Fiber-Epoxy Matrix After Impact 
for 0° Layup Angle 

29 Average Membrane Stress Distribution in the Quarter Plane 
of the Plate for Graphite Fiber-Epoxy Matrix After Impact 
for ±45° Layup Angle 

« 

30 Average Flexural Stress Distribution in the Quarter Plane 
of the Plate for Graphite Fiber-Epoxy Matrix After Impact 
for ±45 Layup Angle, Thickness /Impact Radius Ratio =10.0 

31 Average Flexural Stress Distribution in the Quarter Plane 
of the Plate for Graphite Fiber-Epoxy Matrix After Impact 
for ±45 Layup Angle, Thickness/ Impact Radius Ratio =10.0 

32 Maximum Interlaminar Shear Stress Distribution in the 
Quarter Plane of the Plate for Graphite Fiber- Epoxy 
Matrix After Impact for ±45° Layup Angle 

33 Flexural Mean Stress History, t +t versus Time After 

11 33 

Impact, at Origin for Graphite -Fiber Epoxy 

34 Interlaminar Flexural Stress History A 2 +t 2 

12 32 

Versus Time After Impact for Graphite-Epoxy r/a=0.6 


Note : 


Time/a, (usec/mm) 

for 

Figs. 14-19 

Time/ [b/ (C /p) 1 / 2 )] 
66 

for 

Figs. 20-27 and Figs. 30-32 
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APPENDIX A: THEORETICAL ANALYSIS OF IMPACT IN COMPOSITE PLATES 


Introduction 

In a previous paper, (ref. I) and in section I of this report 
the Author examined the propagation of wave surfaces in composite 
plates, and the response to a line impact load respectively. In 
section II of this report results for the two-dimensional plate 
response to a distributed impact load were presented. This analysis 
makes use of a computational tool called the "Fast Fourier Transform," 
which permits the calculation of the inverse of Fourier transforms on 
the digital computer. The application of this technique to the calculation 
of impact induced stresses is renewed in this section. 

The mathematical model treats the laminated plate as an equivalent 
anisotropic material using a program developed by Chamis (ref. 9) . A modi 
fication of Mindlin's theory of crystal plates is used, which results 
in five two-dimensional stress waves. Two of these waves describe the 
average or membrane stresses, while three other waves are associated 
with the flexural motion. The two former waves are non-dispersive , 
while the flexural waves exhibit strong dispersion at the low frequencies . 
In the Mindlin plate theory, (ref. 7), the impact transverse surface force 
enters the problem directly through the differential equations. To solve 
this initial value problem, a Laplace transform is taken on the time 
\ariable, and a double Fourier transform is taken on the two space 
variables. A solution of this triple transformed problem is obtained 
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in the transform space. Finally, the solution is completed by an 
analytical inversion of the double Fourier transform using the fast 
Fourier transform algorithm. 

The displacement variables used in the theory are described as 

follows: u° and u° represent the midplane displacements in the plane 

of the plate, (see Fig. 9 ); u° represents the transverse midplane 

displacement u 1 and u 1 are a measure of the rotation of a line 
1 3 

normal to the midplane. The equations, which were derived in Reference 1, 
are listed below, where q 2 is a transverse tension force on the plate 
surface and the constants C.^ are the equivalent elastic stiffnesses 
for the composite plate. 
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where 


C 11 = C ll“ C 12 //C 22 


C 33 = C 33- C y C 22 


C 1 3 “ C 13 C 12 C 32 //C 22 


C c , = k C 
66 " 66 

K = 7T*/12 

Analytical Part of Solution: Midplane Motion 

Solutions for the transforms of u° 

l 

easily found. The Laplace variable is denoted by s and the 

Fourier variables are k 1 = k cos a, k 3 = k sin a . The vector formed 


and u° are 
3 
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from (k j , k 3 ) = k represents a Fourier wave number vector corresponding 
to the harmonic frequency a = -is . The resulting expressions become, 

T[L[u 0 ]] = 0 = i q {(A 0 „ + s 2 /k 2 )a_ - A 1? a 0 sin 2 a} cos a/kA 

l 1 22 1 z 

T[L[u 0 ]] = U = i q { (A + s 2 /k 2 )a 0 - A.. a, cos 2 a} sin a/kA (A- 3) 

3 3 11 z 1 

A = (A,. + s 2 /k 2 ) (A„ 0 + s 2 /k 2 ) - A 2 cos 2 a sin 2 a 
1 1 12 

where , 


A }1 = C 1: cos 2 a + C 55 sin 2 a , A 22 = C 33 sin 2 a + C 55 cos 2 a, 


A = c + C 
12 L 13 55 


Cj 2 /2C 22 ’ - ^23^^22 

(A bar indicates a Laplace transform, and q indicates a double Fourier 
transform.) These solutions have the form. 

A 

0 = i q P(s 2 )/kA(s 2 ) 

2 a 2 2 2 

where P(s ), A(s ), are polynomials in s . A (s ) has four zeroes 
in the complex s plane for each k_. These roots have the form, 

s = ± ivj k , ± iv 2 k 
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where v^ and v are the plane wave velocities corresponding to the 

wave normal k . Values of v and v versus a were reported in 
- 12 

Reference 1. The Laplace inversion of (A-3) can then be done by use of 
the convolution theorem; 


U(k,a,t) 


G(t) = 


1 

k 


rt 


q(k,a,t-x) G(x)dT 


(A-4) 


k P(-k 2 v 2 ) sin k v x t k P(-k 2 v 2 ) sin k v t 


Cv^-v 2 ) 


l ' 


(v 2 -V 2 ) 


Analytical Part : Flexural Motion 

By a similar procedure, one can solve for u° , uj , u* , which 
have the form 

W = q (A-5) 

A 1 (s 2 ,k 2 ) 

However in this case the roots of =0 are not proportional to 
k . This means that the velocity of the waves depends on the wavelength. 

It is known from the dispersion relations for these plates (ref. 7) , that for 
k real, there will be six pure imaginary roots of A^s 2 ) = 0; 


s = ±ia> 1 fk) , ±io) 2 Q<), ±iu> 3 (k) 

The Laplace inversion of (A-5) again makes use of the convolution 


theorem and the residue calculus to invert 1/A(s 2 ). Thus we obtain 



(A-6) 


ft 

= q(Jc,t-T) H(r)dT 
1 o 


where 


R(u) 1 ) sin to t R(to ) sin to 2 t 

H(t) = + — 

(u)2-a>2) (to^-to*)^ (ce|-a3|) (aj|-a>|)u> 2 


R(io 3 ) sin to 3 t 


(to^-top (to2-(o|)to 3 


Numerical Inversion: 


The inversion of the Fourier transforms involves integrals of 


the form 


U(x x ,x 3 ) 



UCk 1# k 3 ) e" 1 ^'^dk i dk 3 


(A-7) 


If significant changes in U(x) occur over distances greater than A , ■ 
then the largest wave number of interest will be 


K = 2tt/A 


Thus we may be satisfied with an approximation to UCx) of the form 


U(*) 


4tt ^ 


U(k) e _1 v’*dk dk, 

v 1 3 


-K 


Shifting coordinates, this becomes 


&C{> = 


^(KXj+KXj) 

4t r 2 


2K 

| j ?)(k 1 -K,k 3 -K) e _1 ^'^dk 1 dk 3 
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This latter integral may further be approximated by the following sum. 


% 

U 


e iCfC*,*Kx 3 ) e N £ CI)-K,k (J)-K) 


tt 2 N 2 


I,J=1 


where 

\ * 

or 


kjCI) = 


2 K 
N 


J* Cl-i)' 


V J > -- 2J s 


u 


- e lK( - 1- N^ x l + X 3- ) Z N Z U(I,J) e'i|^- [C I - 1 ) X 1 + ( J ~ 1 ) X 3 ] 

.9 »t 9 -r -r , N 


7T 2 N 2 


I ,J=1 


(A- 8) 


When x and x are continuous variables, a summation of the type 
13 

above must be performed for each point (x^Xg), which makes the cal- 
culation of the sums impractical for a large number of grid points 

(x. ,x_) . However, if x and x take on certain descrete values, there 
1 i 13 . 

is an algorithm (ref. 8), which makes the calculation of these sums feasible for 
a large number of grid points. This algorithm, known as the "fast Fourier 
transform," takes a sampling matrix of the transform, say A(kj,kj), and 
returns the following sampling matrix of the original transformed function. 


T(L,M) 


Z N Z A(k k ) 
J=1 J=1 1 ■ 




CJ-1) (M-l ) , 
N J 


(A-9) 


This operation is known as a descrete finite Fourier transform. To 

put the above expression into this form we choose for the discrete 

values of x. and x the numbers 
1 3 
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CA- 10) 


Xj(L) = j (L-l) , x 3 (M) = j (M-l) 

Numerical Results 

Several checks were made of the fast Fourier computer routine 
(ref. 8) used in this paper. First, known functions were -transformed 
analytically and inverted numerically. These tests revealed that 
one must work with continuous functions if spurious oscillations 
near points of discontinuity are to be avoided. 

Second, known one -dimensional solutions were checked using the 
numerical transform method, such as the impact of a string and the 
impact of a classical beam (see Section I). All these checks 
revealed very close agreement between the numerical transform output 
and the known analytical result. 

t 

The impact pressure distribution used was the following (see 
Figure 13) 

<2 " - P o t 1 - 2 ^ + Cl)") Sin f 

0 

for r<a, (r 2 = x 2 + x 2 ) and t<x 

13 0 

q 2 = 0, for r>a or t>T Q 

The Hertz contact pressure based on static isotropic elasticity has the 
form 


(A-ll) 

f 
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(A-12) 



The principal difference between (A-ll) and (A-12) lies in the infinite 
slope in (A-12) at r=a . The form of the impact pressure (A-ll) was 
chosen to avoid such infinite slope. A distribution with continuous 
first and second derivatives (as the function (A-ll) exhibits) was dictated 
by a desire to have all stresses continuous (i.e. to avoid shocks) and 
thus avoid spurious oscillations in the numerical inversion. This require- 
ment results from the fact that in the Mindlin theory the average mid- 
plane stresses, having as wave sources terms proportional to Vq 2 , have 
the form 




kq 2 (kr)sin k v(t-t)dx^ (a,B 


1,3) 


where v denotes either the first or second wave speed. Thus when 
q^ has the form of a short duration impact i.e. 


q 2 = Q(k) <5 (t) 


the Fourier transforms of the stresses are proportional to 


t n0 'v k 2 Q(k) sin k vt 

ay 
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or for small times. 


9 2 Q(r) 

3r 2 


Thus for the stresses to be continuous at early times after the wave 
arrival, the spacial part of q 2 must have continuous second derivatives, 
which led to the choice of (A- 11) .This conclusion was also reached in numeri- 
cal tests of the fast Fourier program when non-smooth load distributions 
were used. (See Section II of this report). 

Using a "fast Fourier" computer routine written in Fortran IV, 
the induced impact stresses were calculated on both an IBM 360-91, and 
IBM 7094 computers. The grid used was 32 x 32 or 1024 points in the 
x 1 - x a plane. Execution time on the 360-91 was of the order of 12 
seconds, and 22 seconds on the 7094. A denser grid of 64 x 64 was 
also tried with a running time of less than a minute on the IBM 360-91. 

The output data consisted of a matrix (32 x 32) of stress values 
for the quarter plane of the plate. Interpolation, contour plotting, 
and three-dimensional plot routines, developed at the NASA Lewis Research 
Laboratory for use with a "Ca^ Comp" plotter, were used to obtain stress 
contours and three-dimensional plots as shown in Figure 2 of Section II. 

The significant stress levels all lie within the surface bounded 

by the theoretical wave surface. In Figures 14, 17 the 

average or membrane mean stress contours 4-(t° + t° ) for eranhite 

2 ll 33 

fiber/epoxy matrix laminate plates are shown for layup angles of 
0°, ±45°. 
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The stresses shown correspond to the faster wave speed (quasi com- 
pressional wave) which is anisotropic, as is seen from the wave surfaces 
in Figure 10 . The stresses associated with the second wave speed 
(or quasi shear wave) were, much smaller than those in the faster wave 
and are not presented. 

The flexural or bending motion has three waves associated with it. 

The largest stresses however were found in the lowest flexural wave 
which travels at an isotropic wave speed given by 

v 3 = t C 66 K /P] 1/2 

(k = ir 2 /12 , is Mindlin's correction factor (?ef. 3). Stress contours for 

wave are shown in Figures 20- 
25 for graphite fiber/epoxy matrix laminate plates under the 
transverse impact pressure (A-ll). Note that the wave 
front is circular since v 3 is isotropic for laminate plates. Stresses 
in the second and third flexural waves were found to be small. Three- 
dimensional computer plots are shown in Figures 26-32. 

The maximum stress levels were found to occur immediately after 
the end of impact and appeared to propagate along the fiber directions, 
given by the layup angles. 


the mean flexural stress 


1 t l 

2 11 


+ t 


33 


in this 
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APPENDIX B: LIST OF SYMBOLS 


a 



C. . 
ij 

D (I) 

e. . 

E 


F 


k 2 

L[f] 

M 

k 

P 

n 

£ 


t 


t. . 
ij 



o 


u 

u 


2 

1 

1 * 


u; 


half width of impact contact line (one dimensional) or radius 
of impact contact circle (two dimensional) 

ratios of elastic constants 

acoustic tensor components 

half thickness of plate 

anisotropic elastic constants 

column matrix of discrete Fourier Transform 

elastic strain tensor 

Young's Modulus 

impact force 

wave number 

Hertz contact law constant 

Laplace transform. 

mass of impacting object 

wave normal 

Legendre polynomials 

impact loading function 

position vector 

Laplace transform variable 

time 

stress tensor 
displacement vector 

inplane displacements 
transverse plate displacement 
flexural rotation displacements 
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V 


V 


v W X 3 


higher order plate displacements 
wave velocity- 

velocity of impacting object 
cartesian coordinates 


a 

a 


A, A 
<P 


1 


K 

X 

v 


P 

w ,Q 

o 

T 

e 


angle of one- dimensional wave,normal 
relative approach of impacting object and plate 
characteristic acoustic determinant 
fiber layup angle 

Mindlin correction factor, or wave number (see the text) 

wave length 

Poisson's ratio 

density 

frequency 

impact time or time parameter 

one -dimensional wave normal direction 
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APPENDIX C: MATERIAL PROPERTIES 


All the calculations in this report are for the composite material 
consisting of graphite fibers in an epoxy matrix. 

The values of the elastic constants for graphite/epoxy, for various 
ply layup angles, were obtained from an analysis by Chamis (ref. 9). 

The assumed properties of the graphite fibers and epoxy used by Chamis 
are as follows; 

Epoxy Young's Modulus E = 0.57 10 6 psi 

Poisson's Ratio v = 0.36 

Graphite Fiber (Thomel 50) 

"1" Axis along the fiber 

E =50 10 6 psi 
1 1 

E = E =1.0 10 6 psi 

22 33 

v = v =0.20 

12 13 

v = 0.25 

23 

G = 1.3 10 6 
12 

G = 0.7 10 6 
23 

The values of the elastic constants for the composite are given in 
the following table. 
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TABLE I. - STRESS-STRAIN COEFFICIENTS FOR 55 PERCENT GRAPHITE 

FIBER -EPOXY MATRIX COMPOSITE 
[All constants to be multiplied by 10 psi; data obtained from ref . 7. ] 


0° Layup 


27.95 0.3957 0.3957 

1.170 0.4601 

1.170 


0.3552 


±30° Layup 

16.48 0.4118 5.167 0 

1.170 0.4400 0 

3.093 0 

0.3552 


0 

0 

0 

0 

0.7197 


0 

0 

0 

0 

5.491 


0 

0 

0 

0 

0 

0.3552 


0 

0 

0 

0 

0 

0.3552 







FIGURE I • DIAGRAM OF COMPOSITE PLATE 
GEOMETRY AND NOTATION 








WAVE NUMBER kb 

FIGURE 2 -DISPERSION RELATION FOR FLEXURAL MOTION cub vs kb 
FOR 55% GRAPHITE- EPOXY ± 45° LAYUP ANGLE 
FOR WAVE NORMAL a = 0 ALONG X, AXIS 
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x/a 

FIGURE 3 • a) WAVE IN AN INFINITE STRING FOR AN INITIAL TRIANGULAR DISPLACEMENT 
(ONLY RIGHT HALF SHOWN) 

b) COMPUTER CALCULATION OF DISPLACEMENT WAVES IN AN INFINITE STRING 
DUE TO IMPULSE LOADING EQU. (19) 
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FIGURE 5- MIDPLANE DISPLACEMENT WAVES, 55% GRAPHITE- EPOXY COMPOSITE , ± 45° LAYUP 
ANGLE a/b = 10.0 , IMPACT LOAD ALONG X 3 AXIS 
- COMPARISON OF MINDLIN AND CLASSICAL THEORIES 




NORMALIZED DISPLACEMENT UjCg^Pob 



FIGURE 6 • CENTER MIDPLANE DISPLACEMENT vs TIME 55% 

GRAPHITE - EPOXY COMPOSITE , ± 45° LAYUP ANGLE 
IMPACT LOAD ALONG X 3 AXIS , a/b= 10 , COMPARISON 
OF MINDLIN AND CLASSICAL THEORIES 






FIGURE 7- CENTER MIDPLANE DISPLACEMENT vs LAYUP ANGLE 55% 
GRAPHITE FIBER - EPOXY MATRIX COMPOSITE a/b = 10.0 


NORMALIZED MEAN MEMBRANE STRESS 



FIGURE 8- MEAN INPLANE STRESS (tfj+t° 3 )/2 vs DISTANCE x/o 55% GRAPHITE 
FIBER - EPOXY MATRIX COMPOSITE PLATE FOR VARIOUS LAYUP 


ANGLES, LOAD ALONG x 3 AXIS 
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FIGURE 10. a) WAVE SURFACES FOR GRAPHITE /EPOXY 0° and ±15° FIBER 
LAYUP ANGLE 



t = 10- 6 sec. 
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HERTZIAN IMPACT CONTACT TIME FOR GRANITE AND ICE SPHERES 

ON GRAPHITE/ EPOXY COMPOSITE 


FIGURE 11 





MAXIMUM IMPACT PRESSURE I0 3 psi 



figure 12a MAXIMUM HERTZIAN IMPACT PRESSURE 

FOR GRAPHITE /EPOXY COMPOSITE 



PRESSURE, 10 8 NEWT./m 2 




MAXIMUM IMPACT FORCE 10° NEWTONS 



figure 12b MAXIMUM HERTZIAN IMPACT FORCE 
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FIGURE 14 
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FIGURE 15 
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MEMBRANE STRESS { (t,, +133) 
GRAPHITE FIBER/EPOXY MATRIX 



FIGURE 16 
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MEMBRANE STRESS +*33) 
GRAPHITE FIBER/EPOXY MATRIX 

\YUP ANGLE : ± 45° 

DRM. IMPACT TIME T = 0.l t NORM. TIME : 
RELATIVE STRESS LEVELS 












FIGURE 19 









FIGURE 22 









INTERLAMINAR SHEAR STRESS 
GRAPHITE FIBER/ EPOXY MATRIX 
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FIGURE 24 





INTERLAMINAR SHEAR STRESS ■/*(♦, i-M&> 
GRAPHITE FIBER/EPOXY MATRIX 



FIGURE 25 
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TRANSVERSE DISPLACEMENT 
GRAPHITE FIBER/ EPOXY MATRIX 



FIGURE 27 


MEMBRANE STRESS \ (t„ + t 33 ) 
GRAPHITE FIBER/ EPOXY MATRIX 



FIGURE 28 


Mean In-Plane Stress Profile After Impact 




FIGURE 29 


FLEXURAL STRESS 5 (t|, + t 33 ) 
GRAPHITE FIBER/EPOXY MATRIX 
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FIGURE 30 




FLEXURAL STRESS 
GRAPHITE FIBER /EPOXY MATRIX 






STRESS IN UNITS OF Qo/SeelC^psi 



FLEXURAL STRESS t,, + t 33 AT ORIGIN vs TIME, a/b = 1.0 
NORMALIZED CONTACT TIME TV 3 /b = I.O (V 3 =I.I8 mm/^isec) 


FinURK 33 
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NORMALIZED CONTACT TIME rV 3 /b = 1.0 (V 3 =I.I8 mm/^.sec 
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